#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "nrsrc/nr.h"
#include "nrsrc/nrutil.h"

#include "prototypes.h"
#include "globvars.h"
#include "cooling.h"



void set_sfr_parameters(void)
{
  PhysDensThresh = SF_THRESH_MinPhysDens_HpCM3 * PROTONMASS / HYDROGEN_MASSFRAC / UnitDensity_in_cgs;

  PhysEngyThresh = SF_EOSEnergyAtThreshold_ERG * UnitMass_in_g / UnitEnergy_in_cgs;

  PhysPresThresh = GAMMA_MINUS1 * PhysDensThresh * PhysEngyThresh;

  SchmidtLawCoeff = SF_SchmidtLawCoeff_MSUNpYRpKPC2 *
    pow(pow(CM_PER_MPC / 1.e6, 2) / SOLAR_MASS * sqrt(GasFraction), SchmidtLawExponent - 1) / (1.e6 * SEC_PER_YEAR);

  SchmidtLawCoeff *= (pow(UnitLength_in_cm, 2) * UnitTime_in_s / UnitMass_in_g);

  printf("PhysDensThresh = %f, PhysEngyThresh = %f, PhysPresThresh = %f, SchmidtLawCoeff = %f\n",
	 PhysDensThresh, PhysEngyThresh, PhysPresThresh, SchmidtLawCoeff);
}

void init_clouds(void)
{
  set_sfr_parameters();

  U4 = PhysEngyThresh;

  //  PhysDensThresh = SF_THRESH_MinPhysDens_HpCM3 * PROTONMASS / HYDROGEN_MASSFRAC / UnitDensity_in_cgs;

  //  printf("Computed: PhysDensThresh = %g  (int units)         %g h^2 cm^-3\n", PhysDensThresh,
  //	 PhysDensThresh / (PROTONMASS / HYDROGEN_MASSFRAC / UnitDensity_in_cgs));

  dump_eqn_of_state();
  integrate_surfacedensity();
}


void dump_eqn_of_state(void)
{
  double rho;
  FILE *fd;

  fd = fopen("eqs.txt", "w");
  for(rho = PhysDensThresh / 1e4; rho < PhysDensThresh * 1e8; rho *= 1.1)
    {
      fprintf(fd, "%g %g\n", rho / (PROTONMASS / HYDROGEN_MASSFRAC / UnitDensity_in_cgs), eqn_of_state(rho) / (GAMMA_MINUS1 * rho) / UnitMass_in_g * UnitEnergy_in_cgs);
    }
  fclose(fd);
}



double eqn_of_state(double rho)
{
  double P;

  if(rho > PhysDensThresh)
    P = PhysPresThresh * pow(rho / PhysDensThresh, SF_EOSGammaEffective);
  else
    P = GAMMA_MINUS1 * rho * PhysEngyThresh;

  return P;
}


double get_gamma(double rho)
{
  double gamma;

  if(rho > PhysDensThresh)
    gamma = SF_EOSGammaEffective;
  else
    gamma = 1;

  return gamma;
}




void integrate_surfacedensity(void)
{
  double rho0, rho, rho2, q, dz, gam, sigma = 0, sigma_u4, sigmasfr = 0, P1;
  double sigmacheck, rhoold;
  double x = 0, P, P2, drho, dq;
  double meanweight, u4, z, tsfr;
  double rhomin, rhomax;
  int bin, rep;
  FILE *fd;

  if(N_GAS == 0)
    return;

  u4 = PhysEngyThresh;

  NumGasScaleLengths = -log((GAMMA - 1) * u4 * H / G / MaxGasDiskHeight / (GasFraction * M_DISK));

  if(NumGasScaleLengths < 8.0)
    NumGasScaleLengths = 8.0;

  printf("NumGasScaleLengths = %g\n", NumGasScaleLengths);

  printf("integrating surface density...\n");

  fd = fopen("surface_density.txt", "w");


  rhomin = PhysDensThresh / 100;
  rhomax = 10000 * PhysDensThresh;


  for(rep = 0; rep < TABSIZE; rep++)
    rhocentral[rep] = exp(log(rhomin) + rep * (log(rhomax) - log(rhomin)) / (TABSIZE - 1));

  for(rep = 0; rep < TABSIZE; rep++)
    {
      rho0 = rhocentral[rep];

      z = 0;
      rho = rho0;
      q = 0;
      dz = 0.0025;

      sigma = sigmasfr = sigma_u4 = 0;

      while(rho > 0.00001 * rho0)
	{
	  P = eqn_of_state(rho);
	  gam = get_gamma(rho);

/* 	  if(rho > PhysDensThresh) */
/* 	    { */
/* 	      gam = SF_EOSGammaEffective; */
/* 	    } */
/* 	  else */
/* 	    { */
/* 	      gam = 1.0; */
/* 	      sigma_u4 += rho * dz; */
/* 	    } */

	  drho = q;
	  dq = -(gam - 2) / rho * q * q - 4 * PI * G / (gam * P) * rho * rho * rho;

	  sigma += rho * dz;

	  if(rho > PhysDensThresh)
	    sigmasfr += SchmidtLawCoeff * rho * dz *
	      pow(P * GAMMA / GRAVITY, (SchmidtLawExponent - 1) / 2);

	  rho += drho * dz;
	  q += dq * dz;
	}


      sigma *= 2;		/* to include the other side */
      sigmasfr *= 2;
/*       sigma_u4 *= 2; */


      bin = 0;
      z = 0;
      rho = rho0;
      q = 0;
      drho = 0;

      sigmacheck = 0;

      while(bin < ZSIZE && rho > 0.00001 * rho0)
	{
	  dz = list_z[bin + 1] - list_z[bin];

	  rhoold = rho;

	  sigmacheck += 0.5 * dz * rho;

	  rho += 0.5 * dz * drho;

	  P = eqn_of_state(rho);
	  gam = get_gamma(rho);

/* 	  if(rho > PhysDensThresh) */
/* 	    gam = SF_EOSGammaEffective; */
/* 	  else */
/* 	    gam = 1.0; */

	  drho = -4 * PI * G * sigmacheck * rho * rho / P / gam;

	  rho = rhoold + drho * dz;

	  sigmacheck += 0.5 * rho * dz;

	  bin++;
	}

      sigmacheck *= 2;

      fprintf(fd, "%g %g %g %g %g\n", rho0, sigma, sigmasfr, sigma_u4, sigmacheck);

      sigmatab[rep] = sigma;
    }

  fclose(fd);

  printf("done.\n");
}










/* returns the maximum of two double
 */
double dmax(double x, double y)
{
  if(x > y)
    return x;
  else
    return y;
}

/* returns the minimum of two double
 */
double dmin(double x, double y)
{
  if(x < y)
    return x;
  else
    return y;
}



/* returns the maximum of two integers
 */
int imax(int x, int y)
{
  if(x > y)
    return x;
  else
    return y;
}

/* returns the minimum of two integers
 */
int imin(int x, int y)
{
  if(x < y)
    return x;
  else
    return y;
}
